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ABSTRACT 

We present evidence that relativistic shocks propagating in unmagnetized plasmas can self- 
consistently accelerate particles. We use long-term two-dimensional particlc-in-ccll simulations to 
study the well-developed shock structure in unmagnetized pair plasma. The particle spectrum down- 
stream of such a shock consists of two components: a relativistic Maxwellian, with characteristic 
temperature set by the upstream kinetic energy of the flow, and a high-energy tail, extending to 
energies > 100 times that of the thermal peak. This tail is best fitted as a power law in energy with 
index —2.4 ± 0.1, modified by an exponential cutoff. The cutoff moves to higher energies with time 
of the simulation, leaving a larger power law range. The number of particles in the tail is ~ 1% of 
the downstream population, and they carry ~ 10% of the kinetic energy in the downstream. Upon 
investigation of the trajectories of particles in the tail, we find that the energy gains occur as particles 
bounce between the upstream and downstream regions in the magnetic fields generated by the Weibel 
instability. We compare this mechanism to the first order Fermi acceleration, and set a lower limit on 
the efficiency of shock acceleration process. 

Subject headings: acceleration of particles - collisionless shocks - gamma rays: bursts 



1. INTRODUCTION 

Acceleration of particles in collisionless shocks is at the 
heart of most models of nonthermal phenomena in the 
Universe. Observations of synchrotron emission from as- 
trophysical sources suggest that collisionless shocks in 
pulsar wind nebulae (PWNe), jets from active galactic 
nuclei (AGNs), gamma-ray bursts (GRBs) and super- 
nova remnants (SNRs) can convert a significant fraction 
of the flow energy into relativistic particles with power 
law nonthermal spectra. This is usually attributed to 
the first-order Fermi mechanism - a process in which 
particles scatter between the upstream and downstream 
regions of shocks to gain energy (see Blandford & Eich- 
ler 1987 for review). Yet, despite its significance in as- 
trophysics, the Fermi mechanism has not been demon- 
strated to work self-consistently from first principles, and 
its efficiency and conditions for operation are not well 
understood. Most progress in studying shock acceler- 
ation has been made using Monte Carlo test particle 
simulations (e.g., Ostrowski & Bednarz 1998, Ellison & 
Double 2004) and semi-analytic kinetic theory methods 
(e.g., Kirk et al. 2000, Achterberg et al. 2001, Keshet 
& Waxman 2005). Both methods study particle accel- 
eration with certain assumptions about the scattering 
processes near shocks, which in turn depend on the na- 
ture of magnetic turbulence in the flow. Whether re- 
alistic shock turbulence leads to particle acceleration is 
currently unknown. This is especially true in the case 
of relativistic shocks, where observational constraints on 
the turbulence properties are lacking. The injection ef- 
ficiency, or the fraction of particles in the flow that are 
nonthcrmally accelerated, may depend on the details of 
the shock transition region and is not well constrained. 
Resolving these issues requires a self-consistent calcu- 
lation that can simultaneously capture the physics of 
the shock-generated turbulence and the particle accel- 
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eration processes. In this Letter, we demonstrate via ab- 
initio particle-in-cell (PIC) simulations that relativistic 
collisionless shocks propagating in initially unmagnetized 
electron-positron pair plasmas naturally produce accel- 
erated particles as part of the shock evolution. These 
particles form a power law tail, and we argue that the 
acceleration process is very similar to Fermi acceleration. 
We are thus able to measure the efficiency of shock accel- 
eration without any assumptions. The particular case of 
unmagnetized relativistic shocks that we consider here 
is relevant to the prompt and afterglow emission from 
GRBs (e.g., Waxman 2006). This is also one of the clean- 
est tests of the shock acceleration theory because in an 
initially unmagnetized flow all magnetic turbulence has 
to be self-generated. In <|2] we describe the simulations 
and shock structure and then discuss the particle accel- 
eration mechanism in Sj3] 

2. COLLISIONLESS SHOCK STRUCTURE 

Collisionless shocks in unmagnetized relativistic flows 
are mediated by the Weibel instability (Weibel 1959; 
Medvedev & Loeb 1999; Gruzinov & Waxman 1999), 
which converts the free energy of anisotropic stream- 
ing of interpenetrating flows into small scale (skindepth) 
magnetic fields. The fields grow to subequipartition lev- 
els and deflect and randomize the bulk flow, creating 
the shock compression (Kato 2005, Milosavljevic et al. 
2006). PIC simulations of colliding relativistic shells 
have confirmed the general picture of the initial stages 
of the instability (Kazimura et al. 1998; Nishikawa et al. 
2003, 2005; Silva et al. 2003; Frederiksen et al. 2004; 
Hededal et al. 2004), and have been evolved through 
the shock formation stage in both 2D (Gruzinov 2001; 
Kato 2007; Spitkovsky 2008, hereafter S08; Chang et 
al. 2008, herefter CSA08) and 3D (Spitkovsky 2005, 
hereafter S05). A persistent feature of all unmagne- 
tized shock simulations had been the downstream par- 
ticle spectrum that is well approximated by a relativistic 
Maxwellian (S05; Kato 2007). Reports of nonthermal 
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spectral components in short 3D simulations (Nishikawa 
et al. 2003, Hededal et al. 2004) can probably be at- 
tributed to the incomplete flow thcrmalization before a 
shock fully forms. Long simulations, as presented in this 
Letter, show that additional nonthermal spectral compo- 
nents develop over time. 

We use the electromagnetic PIC code TRISTAN-MP 
(Buncman 1993; S05) to simulate a rclativistic shock 
propagating through an unmagnetized e plasma. The 
shock is triggered by reflecting an incoming cold "up- 
stream" flow that propagates with Lorentz factor 70 = 15 
in the —x direction from a conducting wall at x = 0. The 
simulation is performed in the "wall" or "downstream" 
frame (SOS; S08). To maximize feasible simulation size 
we use the 2D version of TRISTAN-MP with flow direc- 
tion lying in the simulation plane. Although we track all 
three components of velocity, in the case of an unmagne- 
tized 2D shock only the in-plane velocities and currents 
are excited. The only magnetic field component is then 
the out-of-plane B z (CSA08). 

Our numerical parameters are as follows: the 
skindepth based on the upstream density n\ is c/u> p — 
(47re 2 ni/7om e c 2 ) -1 / 2 = 10 cells, the transverse box size 
l y = 400c/(jj p , and timestep ui p At — 1/22, with 8 par- 
ticles/cell/species in the incoming flow. We have also 
tried larger box sizes (l v < 1024c/wp) and more parti- 
cles/cell/species (up to 32) with no significant changes 
to the results. The simulations are evolved for times 
> 10 4 w- l . The final box length is l x > 10 4 c/<*V 

In Fig. [T|we show the internal structure of the shock at 
time uj p t = 8400. Density (Fig. Ilk) and magnetic field 
energy (Fig. [Tja) show the Weibel filamentation in the 
upstream region, with characteristic transverse filament 
scale ~ 20c /ui p . The filaments are advected with the 
upstream flow into the shock and undergo merger. The 
flow rapidly decelerates over ~ 50c/cj p , and the density 
increases to nz/ni = 3.13 (Fig. uk), as given by the 
jump conditions for a 2D relativistic gas measured in the 
downstream frame (CSA08; S08). The shock moves to 
the right at 0.47c (shock Lorentz factor with respect to 
the upstream is T shock = V^lo = 21.2). Transversely 
averaged magnetic energy reaches a peak at ~ 15% of 
equipartition and then decays (Fig. [Tji; CSA08). In 
Fig. [l^, we show the longitudinal momentum phase 
space x—jP x for electrons, where the phase space density 
of particles is shown as a 2D histogram (positron phase 
space is nearly identical). The incoming stream appears 
at a negative 4- velocity jfl x « —15, and the flow stops 
and thermalizes after the shock. In the upstream, there 
is a clear population of returning particles with positive 
4- velocity (Milosavljevic et al. 2006; Kato 2007; S08). 
Such hot, low density particles moving in front of the 
shock drive filamentation in the upstream and arc essen- 
tial for maintaining the shock as a self-propagating struc- 
ture far from the wall. Phase space snapshots 7/3 x — jf3 y 
at three cross sections through the flow are shown in 
Fig. [ijj-i. The downstream distribution displays near 
isotropy, while in the upstream the returning particles 
are anisotropic, with an excess of high-energy particles 
that move at large angles to shock normal. 

In Fig. [l]-l we show particle spectra in the down- 
stream and the upstream regions, measured in the ref- 
erence frame of the downstream medium. In the up- 
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Fig. 1. — Internal structure of a relativistic pair shock at uj p t = 
8400: a) Density in the simulation plane (scaled by the upstream 
density); b) Magnetic energy density eg = £? 2 /47r7oni mc 2 ; c) 
Transversely averaged plasma density; d) Transversely averaged 
£b; e) Longitudinal phase space for electrons; f) Transverse phase 
space for electrons; g-i) Momentum space f/3 x — ■y/3 y at three slices 
through the flow (locations marked by arrows); j-1) Electron spectra 
in corresponding slices. 

stream, the returning particles appear as a high-energy 
peak in addition to the incoming flow. In the down- 
stream, the spectrum consists of the main peak, cen- 
tered around 7 *=s 12, and a high energy "tail" extending 
to 7 s» 1000. This tail develops in all our runs at times 
ujpt > 10 3 . In Fig. p] we show a fit to the spectrum 
from a slice centered at 500c/w p downstream from the 
shock at time uj p t = 10 4 . The fit consists of a rela- 
tivistic 2D Maxwellian plus a power law with an expo- 
nential cutoff of the form f(j) = Cijexp(— 7/A71) + 
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C 2 7 _Q min[l,exp(-(7 - 7 ctlt )/A7 cut )], with normaliza- 
tions C\ and C2, such that C2 = for 7 < 7 m j„. In 
the fit, the Maxwellian spread is A71 = 6 (close to but 
smaller than A7 = (70 — l)/2 = 7 expected from the 
full thcrmalization of the upstream flow). The power 
law begins at 7 m i„ = 40 with index a — 2.5, and the 
high energy cutoff starts at j cu t — 300 with a spread 
A7 cu t = 100. For comparison, the red line in Fig. [2]i 
shows a fit with two Maxwellians, demonstrating a clear 
deficit at intermediate energies, while fitting well the low 
and high energy ends. The power law stretch develops 
over time, as can be seen in Fig. [2)3, where a sequence of 
spectra measured at a fixed distance behind the shock is 
shown at different times. The high energy cutoff and 
maximum energy in the tail grow with time and the 
power law range expands. The power law index varies 
between 2.3 — 2.5 (smaller values are obtained if the fit 
is not required to match the transition region between 
the Maxwellian and the power law). At Lo v t — 10 4 , the 
tail at 7 > 75 contains ~ 1% of particles and ~ 10% of 
energy in the downstream. 

3. ACCELERATION MECHANISM 

We studied the mechanism that populates the 
suprathermal tail by tracing the orbits of particles that 
gain the most energy. The main acceleration happens 
near the shock, as seen from the excess of particles with 
large 4-velocity near the shock in Fig. [T^. The space- 
time trajectory x(t) — x s hock{t) and the acceleration his- 
tory 7(f) for four representative particles are shown in 
Fig. [3] The vast majority of particles in the flow go 
through the shock only once and never return to the up- 
stream again (orange line in Fig. [3]). Some, however, 
can cross the shock several times and gain energy. Af- 
ter acceleration near the shock, these particles escape 
into the upstream or the downstream, populating the 
suprathermal tails (red, green and blue particles). The 
particles that gain the most energy (red and blue lines) 
undergo several reflections between the downstream (or 
the shock layer) and the upstream, with the largest en- 
ergy gains coming from reflections in the upstream re- 
gion (Fig. |3j. Upon each reflection these particles gain 
energy AE^ E, as expected in relativistic shocks. In 
Fig. [3] we overplot the transversely averaged magnetic 
energy as line plots stacked in time. Note, that all quan- 
tities are still measured in the downstream frame and are 
only shifted in space so that the shock appears station- 
ary. Magnetic fluctuations associated with the upstream 
filaments carry motional electric field (E y ) as they are 
advected towards the shock. Particles moving against 
the flow in these fields scatter with a net energy gain (in 
contrast, deflections in the downstream region result in 
no energy gain as seen in the downstream frame). The 
shock width of ~ 50c/w p is roughly equal to several Lar- 
mor radii of the thermal particles in the self-generated 
field. The accelerated particles, however, have Larmor 
radii that significantly exceed the width of the shock. A 
high energy particle moving along the shock normal will 
not reflect at the shock, and will leave through the down- 
stream, where the fields are weaker. Instead, the parti- 
cles that get accelerated to the highest energies move 
almost parallel to the shock surface, across the magnetic 
filaments. For them, the alternating magnetic polarity of 
the filaments represents magnetic fluctuations on scales 
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Fig. 2. — Main panel: Particle spectrum in a 100c/u;p-wide slice 
at 500c/oj p downstream from the shock at time uj p t = 10 4 (black 
line with error bars). Red line: a fit with a 2D Maxwellian (yel- 
low dashed) plus a power law (blue dash-dot) with high-energy 
exponential cutoff. Subpanels: a) fit with a sum of high and low 
temperature Maxwellians (red line), showing a deficit at interme- 
diate energies; b) evolution of particle spectrum in a downstream 
slice with time: WOOuip 1 (blue), SSOOuip 1 (green), 10 4 ujp (red). 
Black dashed line shows "-y -2 - 4 power law. 
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Fig. 3. — Left panel: horizontal position as a function of time 
for four representative particles (color lines) overplotted on trans- 
versely averaged profiles of magnetic energy (gray lines). Right 
panel: particle energies (shown with corresponding color lines) as 
a function of time. All horizontal positions are shifted by x s i loc k(t) 
to align them with the shock location. All quantities are measured 
in the downstream frame. 

smaller than the Larmor radius. The deflections in the 
upstream toward the downstream are thus grazing in- 
cidence collisions with moving magnetic islands. The 
deflections towards the upstream happen within several 
hundred skindepths behind the shock, where the mag- 
netic field is strongest during the simulation. 

4. DISCUSSION 

We have shown that relativistic collisionless shocks can 
self-consistently accelerate nonthermal particles. The 
magnetic fields that are created as part of the Weibel 
turbulence near shocks are sufficient to inject particles 
from the thermal pool into the shock acceleration pro- 
cess, and our findings do not require any initial assump- 
tions about the turbulence spectrum. The acceleration 
we observe involves repeated crossings of the shock by a 
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small fraction of particles. In this sense, it is tempting 
to associate this acceleration process with the first order 
Fermi acceleration in shocks. The absence of a coherent 
background field in our problem rules out shock-surfing 
and shock-drift acceleration. Deviations in particle tra- 
jectory are due to interactions with many magnetic fila- 
ments of alternating magnetic polarity and is, therefore, 
diffusive in nature. The exponential cutoff at high en- 
ergy in the downstream spectrum can be attributed to 
the finite acceleration time in the simulation and the fact 
that efficient particle scattering occurs in a layer of finite 
width (Bykov &; Uvarov 1999). We have not reached a 
steady state, however. Both the maximum particle en- 
ergy and the extent of the particle precursor continue to 
grow linearly with time, and the energy in the tail grows 
logarithmically. The region where particles scatter also 
increases. We conclude that the simulations show the 
beginnings of the Fermi acceleration process, and larger 
simulations may display further evolution of the power 
law tail. However, even now we can place a lower limit 
on the efficiency of shock acceleration: at least ~ 1% of 
particles are injected into the nonthermal tail. 

The amount of energy in the tail at the end of our 
simulation is significant enough (~ 10%) to affect the 
jump conditions at the shock. The Maxwellian part of 
the downstream spectrum appears at a lower tempera- 
ture than would be expected from simple thermalization 
of the upstream kinetic energy. Yet, the mean energy 
per particle in the downstream region is still equal to 
the kinetic energy per particle far upstream. The differ- 
ence is due to the contribution from the suprathermal 
tail to the energy balance. In fact, in addition to the 
"cooling" of the Maxwellian component, we observe that 
the fit by a Maxwellian plus power law becomes worse 
with time, mainly because the low-energy peak becomes 
more distorted, with more higher-energy particles than 
predicted by a pure Maxwellian. This transition region 
may require a different spectral component. Extrapo- 
lating this process, we expect that the eventual fraction 
of power law particles may exceed few percent, and the 
thermal peak may get significantly eroded. This process 
may explain the absence of a clear signature of postshock 
thermal emission in PWN and GRB shocks. 

The acceleration process described in this Letter is very 
robust. The development of the tail is not sensitive to 
the transverse size of the box (as long as at least several 
Weibel filaments can fit into the box), or to the vary- 
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ing 7 of the upstream flow (we tried 70 = 2 to 100). 
We have not studied nonrelativistic shocks, because the 
electromagnetic Weibel instability may be less efficient in 
that case, and the shock structure may change. This is 
the subject of future work. In the electron-ion rclativis- 
tic shocks, however, the beginnings of the tail forma- 
tion are already clear in the unmagnetized simulations 
with mass ratio mj/m e = 100, 70 = 15 in S08. It is 
commonly expected that electrons have more difficulty 
being injected into the Fermi process, because their Lar- 
mor radii have to be comparable to the shock thickness 
that is set by ion dynamics. The efficient energy transfer 
between ions and electrons found in S08 eliminates this 
problem, and the relativistic electron-ion shocks thus be- 
have very similarly to electron-positron shocks. There- 
fore, given sufficient computational resources we expect 
to see nonthermal tails in both electrons and protons, 
accelerated by the Fermi process. Although our simula- 
tions are two-dimensional, we expect the same process to 
work in three-dimensional shocks. We were able to run 
3D pair shock simulations for 2000u~ 1 , and have seen 
the beginnings of nonthermal tail formation very similar 
to the 2D case (Spitkovsky & Arons, in preparation). 

Our results underscore the importance of long-term 
PIC simulations for characterizing both the collisionlcss 
shock structure and the nonlinear processes that nat- 
urally appear once a self-consistent shock structure is 
in place. Here we addressed the acceleration of parti- 
cles from the thermal pool. The effect of such acceler- 
ated particles on the evolution of magnetic fields in the 
shock will be presented elsewhere (Keshet et al. 2008). 
The self-consistent picture that emerges from PIC shock 
simulations suggests that the astrophysical scenario of 
shocks depositing ~ 10% of energy into nonthermal par- 
ticles is viable on the microphysical level. In application 
to GRBs our results suggest that prompt and afterglow 
shocks can have nonthermal electron energy fraction of 
e e w 0.1. For PWNe, the mechanism presented here 
can work in shocks near the rotational equator of a pul- 
sar, where, due to the reconnection of oppositely directed 
toroidal fields, the local value of magnetization can be 
small enough to allow Weibel-driven shocks. 
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